

Polyfem Examples¶

Some imports for plotting

In [1]:

import plotly.offline as plotly
import plotly.graph_objs as go

#Necessary for the notebook
plotly.init_notebook_mode(connected=True)

and algebra

In [2]:

import numpy as np

Now we import polyfempy

In [3]:

import polyfempy as pf

Plotting utility¶

In [4]:

def plot(vertices, connectivity, function):
    x = vertices[:,0]
    y = vertices[:,1]

    if vertices.shape[1] == 3:
        z = vertices[:,2]
    else:
        z = np.zeros(x.shape, dtype=x.dtype)

    if connectivity.shape[1] == 3:
        f = connectivity
    else:
        #Convert a tet-mesh into face based triangles.
        f = np.ndarray([len(connectivity)*4, 3], dtype=np.int64)
        for i in range(len(connectivity)):
            f[i*4+0] = np.array([connectivity[i][1], connectivity[i][0], connectivity[i][2]])
            f[i*4+1] = np.array([connectivity[i][0], connectivity[i][1], connectivity[i][3]])
            f[i*4+2] = np.array([connectivity[i][1], connectivity[i][2], connectivity[i][3]])
            f[i*4+3] = np.array([connectivity[i][2], connectivity[i][0], connectivity[i][3]])

    mesh = go.Mesh3d(x=x, y=y, z=z,
                     i=f[:,0], j=f[:,1], k=f[:,2],
                     intensity=function, flatshading=function is not None,
                     colorscale='Viridis',
                     contour=dict(show=True, color='#fff'),
                     hoverinfo="all")
    layout = go.Layout(scene=go.layout.Scene(aspectmode='data'))
    fig = go.Figure(data=[mesh], layout=layout)

    plotly.iplot(fig)

Plate hole example¶¶

Here is the python version of the plate with hole example explained here

Set the mesh path

In [5]:

mesh_path = "plane_hole.obj"

create settings

In [6]:

settings = pf.Settings()

Pick linear $P_1$ elements. If the mesh would be a quad it would be $Q_1$

In [7]:

settings.discr_order = 1

Normalize the mesh to be in [0,1]^2

In [8]:

settings.normalize_mesh = True

Choose Young's modulus and poisson ratio

In [9]:

settings.set_material_params("E", 210000)
settings.set_material_params("nu", 0.3)

We are fine with linear material model

In [10]:

settings.tensor_formulation = pf.TensorFormulations.LinearElasticity

Setup the problem

In [11]:

problem = pf.GenericTensor()

sideset 1 has zero displacement in $x$

In [12]:

problem.add_dirichlet_value(1, [0, 0], [True, False])

sideset 4 has zero displacement in $y$

In [13]:

problem.add_dirichlet_value(4, [0, 0], [False, True])

sideset 3 has a force (Neumann) of [100, 0] applied

In [14]:

problem.add_neumann_value(3, [100, 0])

set the problem

In [15]:

settings.set_problem(problem)

Solve!

In [16]:

solver = pf.Solver()

solver.settings(str(settings))
solver.load_mesh(mesh_path)

solver.solve()

Get the solution

In [17]:

[pts, tets, disp] = solver.get_sampled_solution()

diplace the mesh

In [18]:

vertices = pts + disp

and get the stresses

In [19]:

mises, _ = solver.get_sampled_mises_avg()

finally plot with the above code

In [20]:

plot(vertices, tets, mises)

Torsion¶

In [21]:

mesh_path = "square_beam_h.HYBRID"

settings = pf.Settings()
#It is an hex-mesh so we are using Q1
settings.discr_order = 1
settings.normalize_mesh = False

#Choose Young's modulus and poisson ratio
settings.set_material_params("E", 200)
settings.set_material_params("nu", 0.35)


#Non-linear problem
settings.tensor_formulation = pf.TensorFormulations.NeoHookean
#we want to do 5 steps of incremental loading to avoid ambiguities
settings.nl_solver_rhs_steps = 5

#setup problem with fixed sideset, rotating an number of tours
problem = pf.Torsion()
problem.fixed_boundary = 5
problem.turning_boundary = 6
problem.axis_coordiante = 2
problem.n_turns = 0.5

#coear visualization mesh
settings.vismesh_rel_area = 0.001

settings.set_problem(problem)

#solving!
solver = pf.Solver()

solver.settings(str(settings))
solver.load_mesh(mesh_path)

solver.solve()
#takes ~1 min because it is a complicated non-linear problem!

#get solution and stesses as before
[pts, tets, disp] = solver.get_sampled_solution()
vertices = pts + disp
mises, _ = solver.get_sampled_mises_avg()

#plot the 3d result
plot(vertices, tets, mises)

</div>